PROGRAM test3j

!  Test the RRF program, and specifically the 3j symbol routine, using the
!  sum rules in Brink & Satchler, Appendix II.

!  See data file in tests/test3j.data

!  [ VERBOSE | QUIET ]
!  More or less printing. Default quiet.
!  VERBOSE implies PRINT ALL; QUIET implies PRINT FAILED.

!  PRINT [ ALL | FAILED ]
!  Print all tests, or only those that fail. (None should fail.) Default
!  is to print failed cases only. Printing all results can produce a lot
!  of output.

!  MIN jmin
!  MAX jmax
!  Set the minimum and maximum j values to consider. Default minj=0, maxj=12.

!  SINGLE
!  

!  TEST [[RULE] n]
!  Select and run the test. Only test 1 programmed at present.

USE rrf_module
USE wigner
USE input, ONLY: item, nitems, read_line, reread, readu, readi
IMPLICIT NONE

INTERFACE
  SUBROUTINE start(quiet,primes)
  LOGICAL, INTENT(IN), OPTIONAL :: quiet
  INTEGER, INTENT(IN), OPTIONAL :: primes
  END SUBROUTINE start
END INTERFACE


INTEGER :: a, b, c, d, e, f, p, x, minj=0, maxj=12, test
! INTEGER, PARAMETER :: nprimes=200000
LOGICAL :: eof, verbose=.false., ok, allok, printall=.false., single=.false.
CHARACTER(LEN=20) :: key

write (0,"(a)") "Please wait ..."

call init_rrf(.true.,10000)
print "(a,i0)", "Prime table includes primes up to ", prime(nprime)


do
  call read_line(eof)
  if (eof) exit
  call readu(key)
  select case(key)
    case("VERBOSE")
      verbose=.true.
      printall=.true.
    case("QUIET")
      verbose=.false.
      printall=.false.
    case ("PRINT")
      call readu(key)
      select case(key)
      case("ALL")
        printall=.true.
      case("FAILED")
        printall=.false.
      end select
    case("MIN")
      call readi(minj)
      minj=2*minj
    case("MAX")
      call readi(maxj)
      maxj=2*maxj

    case("SUM-RULE","SUMRULE","SUM", "TEST")
      test=1
      do while (item<nitems)
        call readu(key)
        select case(key)
        case("RULE")
          call readi(test)
        case("SINGLE")
          single=.true.
          exit
        case default
          call reread(-1)
          call readi(test)
        end select
      end do
      select case(test)

      case(1)
        !  Sum rule 1
        !  sum_{ma mb mc md me mf} ( a  b  e )( a  c  f )(  d  b  f )( d  c  e )
        !                          (ma mb me )(ma mc-mf )(-md mb mf )(md-mc me )
        !
        !        x (-)^{d+c+f+md+mc+mf} = { a b e }
        !                                 { d c f }
        !  See my angular momentum lecture notes or Brink & Satchler p. 142
 
        x=2
        print "(/a)", "Sum rule 1"
        if (printall) print "(/a)", "  a     b     e     d     c     f"
        allok=.true.
        !  Arguments a-f are each twice the angular momentum value.
        !  x=2 always but is present for consistency with other routines.
        if (single) then
          call readj(a)
          call readj(b)
          call readj(e)
          call readj(d)
          call readj(c)
          call readj(f)
          printall=.true.
          call sumrule1(a,b,e,d,c,f,x,ok,.false.)
        else
          do a=minj,maxj
            do b=a,maxj
              do e=abs(a-b),min(a+b,maxj),2
                ! args=""
                ! p=1
                ! call put(a, args, p)
                ! call put(b, args, p)
                ! call put(e, args, p)
                ! print "(a)", trim(args)
                ! flush(6)
                do d=minj,maxj
                  do c=abs(d-e),min(d+e,maxj),2
                    do f=abs(a-c),min(a+c,maxj),2
                      call sumrule1(a,b,e,d,c,f,x,ok,.false.)
                      if (.not. ok) then
                        call sumrule1(a,b,e,d,c,f,x,ok,.true.)
                        stop
                      end if
                    end do
                  end do
                end do
              end do
            end do
          end do
        end if
        if (allok) print "(a)", "Test completed -- no errors"

    end select
  end select
end do

CONTAINS

SUBROUTINE sumrule1(a,b,e,d,c,f,x,ok,debug)

!  Arguments a-f are each twice the angular momentum value.
!  x=2 always but is present for consistency with other routines.
INTEGER, INTENT(IN) :: a, b, c, d, e, f, x
LOGICAL, INTENT(OUT) :: ok
LOGICAL, INTENT(IN) :: debug
TYPE(rrf), POINTER :: k=>null(), k1=>null(), k2=>null(), k3=>null(),    &
    k4=>null(), z=>null(), sum=>null()
INTEGER :: ma, mb, mc, md, me, mf
CHARACTER(LEN=35) :: args
CHARACTER(LEN=8) :: verdict
CHARACTER(LEN=200) :: string=""


call setzero(sum)
do ma=-a,a,2
  do mb=-b,b,2
    me=-ma-mb
    if (me<-e .or. me>e) cycle
    do md=-d,d,2
      mc=md+me
      if (mc<-c .or. mc>c) cycle
      mf=md-mb
      if (mf<-f .or. mf>f) cycle
      call threej(a,b,e, ma,mb,me, x, k1)
      call threej(d,c,e, md,-mc,me, x, k2)
      call threej(d,b,f, -md,mb,mf, x, k3)
      call threej(a,c,f, ma,mc,-mf, x, k4)
      call unit(k)
      if (debug) then
        call show3j(a,b,e, ma,mb,me, k1)
        call show3j(d,c,e, md,-mc,me, k2)
        call show3j(d,b,f, -md,mb,mf, k3)
        call show3j(a,c,f, ma,mc,-mf, k4)
      end if
      call mult(k,k1)
      call mult(k,k2)
      call mult(k,k3)
      call mult(k,k4)
      if (modulo((c+d+f+mc+md+mf)/2,2) .ne. 0) call chsign(k)
      if (debug) then
        p=1
        string=""
        call rrf_to_char(k, string, p, 0)
        print "(a,a)", "term:   ", trim(string)
      end if
      call add(sum,k,debug)
      if (debug) then
        p=1
        string=""
        call rrf_to_char(sum, string, p, 0)
        print "(a,a)", "sum:    ", trim(string)
      end if
        
    end do
  end do
end do
call sixj(a, b, e, d, c, f, x, z)
if (debug) then
  call display(sum,1,"sum:")
  call display(z,1,"6j:")
  flush(6)
end if
p=1
args=""
call put(a, args, p)
call put(b, args, p)
call put(e, args, p)
call put(d, args, p)
call put(c, args, p)
call put(f, args, p)
if (equal(sum,z)) then
  verdict="  OK"
  ok=.true.
else
  verdict=" Wrong!"
  ok=.false.
  if (allok .and. .not. printall) print "(/a)", " a    b    e    d    c    f"
  allok=.false.
end if
if (printall .or. .not. ok) then
  p=1
  string=""
  call char4i(sum, string, p, 200)
  p=max(40,p)
  string(p:)=" 6j = "
  p=p+6
  call char4i(z, string, p, 200)
  print "(4a)", args, verdict, "  Sum = ", trim(string)
end if

END SUBROUTINE sumrule1

SUBROUTINE readj(j)
!  Read an angular momentum value m or n/2. The value returned is
!  2m or n respectively

USE input, ONLY: reada
INTEGER, INTENT(OUT) :: j
INTEGER :: p=1, x
CHARACTER(8) :: buffer

call reada(buffer)
p=1
x=2
call number(buffer,p,j,x)
if (x == 1) j=2*j

END SUBROUTINE readj

SUBROUTINE put(j, s,p)

CHARACTER(LEN=*) :: s
INTEGER, INTENT(IN) :: j
INTEGER, INTENT(INOUT) :: p

if (modulo(j,2) == 0) then
  write (unit=s(p:p+2),fmt="(i3)") j/2
else
  write (unit=s(p:p+4),fmt="(i3,a2)") j,"/2"
endif
p=p+6

END SUBROUTINE put

SUBROUTINE show3j(a,b,c,ma,mb,mc, k)

INTEGER, INTENT(IN) :: a,b,c,ma,mb,mc
TYPE(rrf), POINTER :: k
CHARACTER(LEN=80) :: string=""

p=1
string=""
call put(a, string,p)
call put(b, string,p)
call put(c, string,p)
call put(ma, string,p)
call put(mb, string,p)
call put(mc, string,p)
call char4i(k, string, p, 80)
print "(a)", trim(string)

END SUBROUTINE show3j

END PROGRAM test3j
